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On some symmetric multidimensional 
continued fraction algorithms 

Pierre Arnoux * Sebastien Labbe ^ 


Abstract 

We compute explicitly the density of the invariant measure for the Reverse algorithm 
which is absolutely continuous with respect to Lebesgue measure, using a method proposed 
by Arnoux and Nogueira. We also apply the same method on the unsorted version of Brun 
algorithm and Cassaigne algorithm. We illustrate some experimentations on the domain of 
the natural extension of those algorithms. For some other algorithms, which are known to 
have a unique invariant measure absolutely continuous with respect to Lebesgue measure, the 
invariant domain found by this method seems to have a fractal boundary, and it is unclear 
that it is of positive measure. 

Keywords'. Multidimensional continued fractions algorithms. Invariant measure. Natural 
extension. Fractal. 
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1 Introduction 


Many continued fraction algorithms, in one or several dimensions, have a unique ergodic measure 
which is absolutely continuous with respect to Lebesgue measure (the so-called Gauss measure); 
this measure plays an important role in the dynamic of the algorithm, and it is interesting, when 
possible, to give an explicit closed formula for the invariant density. 

It is easy to check if a given function 6 is an invariant density for a map F, since it must be 
solution of the functional equation 5{x) = Y.y,F{y)=x where Jp is the Jacobian of F; but it 

is generally difficult to hnd an explicit solution to this equation. One method, as proposed by 
Arnoux and Nogueira AN93 , is to build a geometric model of the natural extension of F, since 
it is often easier to hnd invariant densities for invertible maps. Another advantage of this method 
is that it gives a suspension how on a manifold which can sometimes be linked to other areas of 
mathematics; for example, for Farey continued fraction, this how can be recognized as the geodesic 
how on the modular surface; for Rauzy induction on interval exchange maps, it is a Teichmiiller 
how. 

The basic idea is as follows. A continued fraction algorithm can often be presented as a 
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piecewise-linear map F on the positive cone of ]R“ (or a subcone of it), given locally as x i— >■ M 
where M is a positive matrix. We can try to build a model F for the natural extension as a skew 
product on a cone of x (x, a) i—)■ (M“^.x, M~^.a). We can look for a cone on which this map 


*Equipe Dynamique, Arithmetique et Combinatoire, Institut de Mathematique de Luminy, CNRS UPR 9016, 
163, avenue de Luminy, Case 907, 13288 MARSEILLE cedex 09, FRANCE, pierre@pierrearnoux.fr 
iRat. B37 Institut de Mathematiques, Grande Traverse 12, 4000 Liege, Belgium, slabbe@ulg.ac.be 


1 





is one-to-one, up to a set of measure 0; contraction arguments (A generalization of Hutchinson’s 
theorem, see AS15|) ensure, in most interesting cases, that such a set is unique, and must exist 


(but it could be of measure 0). 

If such a set D can be found, then F : D ^ D is hj construction a bijection which preserves 
Lebesgue measure and the scalar product x.a. Furthermore, the flow (pt dehned on D by <p 4 (x, a) = 
(e*x, e“*a) also preserves the measure and the scalar product, and commutes with F. Let Di be 
the subset of D dehned by x.a = 1, and be the quotient Di/F; the how ipt projects to the 
quotient, and we can construct a natural extension of the projective map / associated with F as 
hrst return map of the how to a section. 

In this paper, we use this method to hnd the invariant density for two multidimensional con¬ 
tinued fraction algorithms, the Reverse algorithm and the Brun algorithm, and we explain the 
problems encountered to apply this method to some other algorithms. Finding the set D is the 
crucial part of the method. A set D is known for Brun algorithm and is not known for the Hurwitz 
nor the Jacobi-Perron algorithm. Experimentations seem to show that it has a fractal structure 
for Jacobi-Perron and Arnoux-Rauzy-Poincare algorithms. 

In Section as an example, we show how to compute the invariant density for the unsorted 
Farey map. In Sectionwe give the general argument. In Section]^ we apply it to the Reverse 
algorithm, in Section to the Cassaigne algorithm and in Section to the Brun algorithm. In 
Section we explain how numerical computations can show in some cases the invariant set, but 
give more complicated results (fractal sets) in other cases. In Section we make some remarks 
about sorted and unsorted algorithms, acceleration and choices of coordinates, to explain why 
one can hnd many variants of the continued fraction algorithms, and we illustrate that on the 
classical continued fraction. We also give the explicit computation of the invariant density for 
Brun’s algorithm in any dimension. 


2 The unsorted Farey map 

Consider the positive cone A = {(x, y) G | 0 < x, 0 < y} and dehne 


F : {x,y) ^ 


{x, y-x) 
{x-y,y) 


a X < y^ 
ii X > y. 


Remark 1. The map F is well-defined, except on a set of Lebesgue measure 0. There is no 
canonical way to define it on the diagonal x = y; we could choose a convention, but since we will 
be interested in the ergodic properties of F, this is irrelevant. This will occur in all the examples 
of the paper, and we will no more remark about it. 

Since the map F is homogeneous of degree 1 (it commutes with positive homotheties), there is 
an associated projective map /, which can be dehned on the unit simplex. This projective version 
of F on A = {(x, I/) G A I X -|- 1/ = 1} is dehned as / : x i—>■ If we take the hrst variable x 

as coordinate on the unit simplex, / can be dehned as (see Figure 
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Figure 1: The graph of the unsorted Farey map. 


Remark 2. Except for a set of measure 0 (the points {x,y) such that ^ is rational), the positive 
orbit {F^{x,y)} is defined for all n eN, and the norm \x + y\ is strictly decreasing along the orbit; 
the same proof shows that the positive orbit {/”(x)} is well defined for all n E N if and only if x 
is irrational. 


Below we compute an invariant measure of /, using a method that also works for some multi¬ 
dimensional continued fractions algorithms. Note that neither F nor / are bijection. For instance, 
we have F(2, 5) = F(5, 3) = (2, 3): 



One way to transform F into a bijection is to stack the removed part and pile up the result in the 


following way, as done in AN93 
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We thus constructed the following function 


F : {x,y,a,(3) ^ 


{x,y - x,a + (3,(3) 
[x - y,y,a,a + (3) 


ii X < y, 
ii X > y. 


which is an explicit model of the natural extension of F. It is a bijection on the domain D = 
{{x,y,a, (3) e X | x > 0, r/ > 0, a > 0,/3 > 0}; indeed, if we dehne Ai = {{x,y) G | 0 < 
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X <y} and A 2 = {(x, y) G | 0 < y < x}, we have F(Ai x A) = A x A 2 and F(A 2 x A) = A x Ai; 
Since Ai, A 2 is a measnrable partition of A, F sends a partition of D to another. 

Both branches of F can be written linearly as the action by a 4 x 4 matrix of the form 


f M-i 0 \ 
0 ) 


where M 



ii X < y and M 



ii X > y. 


Hence F preserves Lebesgne measnre, since its Jacobian is 1. It also preserves the form xa + yl3, 
that is, the total area of the rectangles. 

We want to hnd a natnral extension and an invariant measnre for the projective map. For this, 
we introdnce the flow (py. 

(ft : {x,y,a,(3) ^ {xe\ye\ae~\ I3e~^). 


This flow has Jacobian 1, thns preserves Lebesgne measnre, and also preserves the form xa + yjS. 
It commutes with F, that is, (pt(F(x)) = F{p>t{x)). 

Dehne now a subset Di of codimension 1 of F by Fi = {(x, y,a, /3) E D \ xa + y/3 = 1}, and 
a surface T, C Di hj T, = {{x,y,a, /3) E Di \ x + y = 1}. 

Di is invariant by F and ipp, F acts without periodic points since it decreases strictly x + y. The 
surface S is a section for ipt, since, for any (x, y, a, (3) E Fi, there exists a unique r = — log(x + y) 
such that ip.r{x,y,a, (3) E E. 

Since (pt commutes with F, we can consider the induced flow on Di/F and the hrst return 
map of this flow to the projection S of the section S on Di/F. By construction, the effect of this 
map on the hrst coordinates is the effect of the map F followed by a renormalization, that is, the 
projective map /. 

We consider the change of coordinates {x,y,a,(3) h-)■ (x,a — /3, —log(x + y),xa + yl3) = 
(x, a', r, e); computation shows that the Jacobian is 1, so Lebesgne measure is given in these coor¬ 
dinates as dx da' dr de. The domain Di is given by e = 1, and the surface of section S by r = 0, 
hence the hrst return map / must leave invariant the measure dxda'. Indeed, a straightforward 
computation shows that / is given in the coordinates of S by 

/: S ^ 

(x, a') I—>■ 

and this map has Jacobian 1. 

We can describe more precisely the surface S by (see Figure]^: 

S = {(x, a', r, e) I r = 0, e = 1, a' -t- /5 > 0, /3 > 0, xa jd = 1} 

= {(x, a', 0,1) I a' -|- 1 — xa' > 0,1 — xa' > 0} 

= {(x, a', 0,1) I - <«'<—} 

X — 1 X 


i-^,1 - X + a'{l - xY^ if X < ^, 
^2 — a'x‘^ — x) if X > ^. 


Therefore the invariant density of / is 


5(x) 



1 

x(l — x) 


One easily check that this density is invariant; a straightforward computation shows that, if 
yi = and y 2 = are the preimages of the point x, we have (5(x) = + ppp)- 
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Figure 2: The domain of the natural extension of the unsorted Farey map. 


Remark 3. That density is unbounded; moreover, the total mass of the measure is infinite. The 
reason is that we have two indifferent fixed points at both extremities of the interval, because we 
have used the additive algorithm. We could accelerate the map around the fixed points by using the 
multiplicative algorithm, and recover the classical finite Gauss measure (more exactly, a variant of 
it, since we use the unsorted map). 

Remark 4. There is one more structure: the map F and the flow ipt leave invariant the symplectic 
form dx A da + dy A dfi, and the flow ipt is the hamiltonian flow associated with the hamiltonian 
xa + yfi. In fact, the set D/F can be identified to the tangent bundle of the modular surface, given 
in coordinates, and (ft is the geodesic flow. 


3 Constructing the measure from the natural extension 


One can extend the method of the previous section to a large range of continued fraction algorithm. 


We recall the method proposed by Arnoux and Nogueira AN93 to hnd heuristics for a geometric 


model of the natural extension of a multidimensional continued fraction algorithm, and an explicit 
formula for the invariant density (Gauss measure) in some cases. 

Let A be a subcone of the positive cone (Two interesting cases are when A = and 
when A = {{xi,X 2 , ■ ■ ■ ,Xd) G | 0 < xi < 0:2 < • • • < Xd} consists of sorted entries). A 
Multidimensional Continued Fraction (MCF) algorithm is a function 


F: A ^ A 

X H- M(x)“Gx. 

where M is an homogeneous function of degree 0, piecewise constant on subcones, that associates to 
each X G A an invertible matrix M(x). Classical references for MCF algorithms are 
In the cases we consider, the entries of M(x) are nonnegative integers. 

To the function F, we can associate its projective version /, dehned in a canonical way from 
the piecewise linear map on the projective space P{A). For explicit computation, one should take 
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coordinates, and a nice way to do that is to consider the codimension 1 compact domain in A, 
A = {x G A I ||x|| = 1} for some norm || ■ ||; it is explicitly given by: 


/: 


A ^ A 


X HG 


iinx)ir 


Assumption 1. We will always suppose that, out of a set of measure 0, we have ||F(a;)|| < ||a; 
this is the case in all our examples. 


Even if there is no canonical coordinates on -P(A), and hence no canonical Lebesgue measure, 
there is a natural Lebesgue measure class, since all changes of coordinates are piecewise projective. 
In many cases, there is a unique invariant (up to a constant) measure in this class, the Gauss 
measure, given by an invariant density. We are interested by explicit formulas for this invariant 
density, and for this, we try to build a geometric model for the natural extension of / with good 
properties. 

A potential geometric model for the natural extension F of the piecewise linear function F can 
be dehned as: 

F \ A X —)■ A X 

M(x)-i 0 V X \ 

0 M{^y y y a y 



To normalize the vector x after the application of F, we dehne the flow ipt as 


ipt'. Ax —)■ Ax 



for all t G M. Note that the following properties are verihed: 

• (pt and F are well dehned, since the matrix M{x) is nonnegative and the dehnition set is a 
product of cones; 

• pt and F both have Jacobian 1; 

• pt and F both preserve the scalar product (x, a); 

• Pt commutes with F\ pt o F = F o pt. 

Remark 5. The map F has no reason to be a bijection, and the main problem is to find a conical 
domain D C A x such that F is a bijection on D. 

From now on, we suppose that we have found such a domain D of positive measure. Then 
Pt and F preserve Lebesgue measure on D, since they are bijections with Jacobian 1. But this 
measure is not very useful, because not only it has inhnite mass, but its projection to A has inhnite 
density. We want to restrict the domain. 

Definition 6. We define Di = {(x, a) G D|(x, a) = 1} 

We can dehne a natural Lebesgue measure on Di by considering coordinates on Di comple¬ 
mented by (x, a) to get a global coordinate system on D, and disintegrating Lebesgue measure. 
The functions F and pt are well-dehned on Di and both preserve locally Lebesgue measure on Di. 
Since they are bijective, they also preserve it globally. 
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Figure 3: The function F applied iteratively on the section S C Di. The gray region corresponds to a 
fundamental domain of the action of F on Di and can be identified to O. 


Definition 7. We define Q = Di/F, the set of orbits of the map F. 

Since ipt and F commute, we can make (ft act on the quotient space Q. We can assimilate Q 
to a fundamental domain of the action of F on Di (this is well defined, since, except for the fixed 
points of F which have measure 0, F decreases the norm of x). 

Definition 8. We define E = {(x, a) G Di \ ||x|| = 1}, and denote by S its projection to Vl. 

Remark 9. There is a canonical projection vr : Di —)■ Di/F, mapping each point on its orbit under 
F. The restriction tt : S —)■ S fs a bijection, up to a set of measure 0: it is onto since we have 
S = 7r(S) by definition, and if it is not one-to-one, there are two distinct elements (x,a), (y, b) G E 
and an integer n > 0 such that (y, b) = F’^(x, a); this implies that ||y|| = ||F’^(x)|| = ||x||, which 
can only happen on a set of measure 0 by assumption. The dynamics happen in the guotient set 
ff, and we are mainly interested in the first return map to E, but all explicit calculus must be done 
in Di and E instead of fl and E. 

Remark 10. The surface E is a section for the flow ipt: for any (x,a) G Di, there is a unigue 
tx = ~log(||x||), depending only on x, such that (pt^(x, a) G E. 

We can define a fundamental domain for the action of F as 

{(y9t(x,a) I (x,a) G E,0 < f < fF(x)} 

(see Figure We consider the first return map / : E —)■ E of the flow on fl. This map can be 
seen in the fundamental domain as the composition of F with the flow: /(x, a) = o -^(x, a). 

Since F preserves Lebesgue measure on Di, there is a well-defined Lebesgue measure on fl, 
which is preserved by ipp, since E is a section, we can define a transverse invariant measure /i on 
E in the usual way: if m is Lebesgue measure on fl and 1/ C E is a measurable set, we define 

_ m{^.(u)\neU.0<s<t} 

t 

From the definition, we have /(x,a) = (/(x), M(x)''^a ■ ||F(x)||), and / can be projected on 
the application /. Hence we can find an invariant measure for / equivalent to Lebesgue measure 
by partial integration on the section E. 
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Proposition 11. AN93 If the measure of T, is positive, by choosing a coordinate system 


x,a 


in which the invariant measure is written dx. da, the invariant measure of f can be computed by 
the formula: 


i(x) = f 


1 da 


{a:(x,a)eS} 


We now describe an explicit change of variables that is useful to compute explicitly all these 
measures in the case of the norm ||x|| = Dehne j/j = Xi for i < d, and r = — log(X]f=i 

dehne bt = Oi — Od ioi i < d and e = ciiXi. A straightforward computation shows that the 
change of coordinates (x, a) i—)■ (y,b,r, e) has Jacobian 1. The domain Di is dehned by e = 1, 
so the invariant measure on 12 is given by nf=i dyt nf=i^ dbi dr; the section E is dehned by r = 1 ; 
since the coordinate r is, up to a sign, the time coordinate of the how, the invariant transverse 
measure for the how, that is, the invariant measure for /, is given explicitly by nf=i dyi n^=i dbi. 


Remark 12. We have made arbitrary choices. The first one was the choice of a coordinates for 
the projective space -P(A); we could choose another section, for example Xd = I (we could also take 
another norm, but it seems in general more convenient to choose an affine section). The second 
one, which is related, is the choice of the section S. The flow (ft and the space 12 are more intrinsic; 
different choices lead to different forms of the same continued fraction; this explains why we can 
find a large range of formulas for apparently similar continued fractions, based on a different choice 
of coordinates. 


4 The Reverse algorithm 


The Arnoux-Rauzy algorithm is a partial algorithm, dehned on the positive cone by subtracting 
the two smallest coordinates from the largest one. It is only dehned if the largest coordinate 
is larger than the sum of the two smaller coordinates, that is, on the complement of the set of 
coordinates which satisfy the triangular inequality. The set of points with inhnite orbits, known as 
the Rauzy gasket AS13 , has measure 0, which make it unusable to hnd rational approximations 


of most points in the positive cone. 

It can be completed in various ways by dehning it on this missing set; for example, one can 
consider the simplest map which sends the central cone onto the positive cone, {x,y,z) h-)■ {—x + 
y + z,x — y + z,x + y — z). 

More formally, we can use the notations of the previous section, and dehne Reverse algorithm 
on the following partition of A = {(x, y, G | a; > 0, ?/ > 0, 2 ; > 0} up to a set of measure zero: 


Ai = {(x, y,z) e A \ X > y + z}, 

A 2 = {(a:, y,z) e A\y > x + z}, 

As = {(a:, y,z) e A \ z > x + y}, 

A 4 = {(x, y,z)EA\x<y + z,y<x + z,z<x + y} 


The name Reverse comes from the fact that its ehect, as a projective map acting on the central 
triangle A 4 of the unit simplex, is an homothety of ratio 2 , followed by a central symmetry which 
reverses the vector with respect to the center of the simplex. 

We dehne the four matrices: 


1 1 1 \ 
0 10 
0 0 1 / 
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Figure 4: Left: the trace on the unit simplex of the partition of A into Aj for i G {1,2, 3,4}. Right: the 
partition of A* into A* for i G {1,2,3,4}. 

and the matrix function M : A —)■ GL{3, Q) such that Mix.) = Ai if and only x G Aj. Recall 
that this matrix function defines the functions F{x) = M(x)”^x, /(x) = F(x)/||F(x)||i and 
F(x, a) = (M(x)“^x, M(x)'''a). 

We show in the left of £gure|^the trace on the unit simplex of the partition of A; each branch 
of the map F sends Aj to all of A. 

Explicitly, the function F is given by 


F : {x,y,z,a,/3,'y) ^ < 


{x — y — z,y, z, a, a + 13, a + 'j) 
{x,y-x-z,z,a + /3,(3,(3 + '^) 
{x,y, z - X - y, a + 'j, 13 + 


^ {—X + y + z,x — y + z,x + y — z, 


/3+7 Q+7 a+y \ 

2 ’ 2 ’ 2 1 


if {x,y,z) G Ai, 
if ix,y,z) G Aa, 
if ix,y,z) G Ag, 
if {x,y,z) G A 4 . 


One can easily find a cone A* such that F is a bijection on F = A x A*: numerical experi¬ 
mentations show that after few iterations of F, (a,/?, 7 ) belongs to the subset of A4 C A of triples 
which satisfy the triangular inequality. Indeed, for any {a,(3,'y) G A, satisfies the 

triangular inequality; and if {a,(3,'j) G A4, then [a,a + (3,a + 7 ) also satishes the triangular 
inequality. Since almost all orbits of F cross A4, with probability 1, after a finite time, the second 
part of the coordinates of F”(x, a) enters A 4 and never escapes. 

Hence we define A* = A4, and a partition of A* by 

a; = {(a, , 3 , 7 ) 6 A-I a < 

a; = {(a,/3,7) 6 A-I/3 < 

A* U R ^ ^ A* I + + 

^3 = e A | 7 <- --}, 

a; = {(a,,3, 7 ) 6 a* I a,/ 3,7 > ^±i±2}. 

Lemma 13. {A{, A2, A3, A4} is a partition of A*, up to a set of measure 0. 
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Proof. It is clear from the definition that A 4 is disjoint from the other 3, since it is defined by 
opposite ineqnalities. Consider a point in the intersection of Aj and A^; we have a,{3 < ^±|± 2 ^ 
hence 7 > which is incompatible with the triangnlar ineqnality. The same proof is valid 

for the two other intersections, hence these fonr sets are disjoint. 

The eqnality a = defines a set of codimension 1 and measnre 0, and similarly for the 

two other ones. For any other point in A*, either all of a, (5 ,7 are larger than which defines 

A 4 , or one (and only one, by the previons paragraph) is smaller, which defines one of the other 
simplices; hence they form a partition, np to a set of measnre 0. We show in the right of fignre 
the trace on the nnit simplex of this partition. □ 

As an immediate conseqnence, we obtain a partition of D\ 

Corollary 14. Up to a set of measure 0, we can write the set D = A x A* as a disjoint union 

4 4 

D = JjAiX A* = JjAx A* 

i=l i=l 

Hence, to prove that F is a bijection on D, since each branch of F is a non-degenerate linear 
map, it snffices to prove: 

Lemma 15. We have F(Aj x A*) = A x A*. 

Proof. We check directly by compntation that F(Aj) = A. We want to show that if x G Aj, then 
Aja. G A*. We nse the notation {x', y', z', a', (3', ■j') = F{x,y, z,a, (3,'y) and we proceed case by 
case. 

If X G Ai, then a' + (3' + = 3a + (3 + 'y > Aa = 4a'. Therefore F(x, y, z, a, /3, 7 ) G A x Aj. A 

similar proof applies for A2 and A3. 

If X G A4, then a' + /3' + = a + f3 + j < 2a + 2(3 = 47 '; a similar proof applies to a' and (3', 

which proves that F{x, y, z, a, (3,^) G A x Aj. 

We have proved inclusion; to prove equality, it is enough to show that the matrix Aj sends the 
extreme points of A* to the extreme points of A*, which is done by computation. See Figure□ 



Figure 5: The function F : Aj x A* —A x Aj is one-to-one and onto for each i G {1,2,3,4}. 


We now apply the results of the previous section. Recall that 

Di = {(a:, y, z, a,(3,'y) G F \ xa + y(3 + Z'y = 1}. 
and the surface of section is defined by 

S = {(x, y, z, a, (3, 7 ) G \xa + yf3 + Z'y = l,x + y + z = l,a + (3>'y,a + 'y>(3,(3 + 'y> a}. 
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It is convenient to take the coordinates on Di defined in the previous section: We keep variables 
X and y and change the coordinates a' = a — ^,f3' = f3 — 'y^e = xa + yf3 + z'y^T = — log(a; + y + z); 
it is readily checked that the Jacobian is one. Since the domain Di is defined by e = 0, and 
coordinate r is the return time of the flow to S, we get 

S = {(x, y, a, /S', r, e) G | r = 0, e = 1, a' + /S' + 7 > 0, a' + 7 > /S', /S' + 7 > a'} 

Furthermore, the invariant measure for the return map / of the flow to S is dx dy da' d/3'. From 
this, we obtain: 

Proposition 16. The density function of the invariant measure 0 / / : A —^ A for the Reverse 
algorithm is 

1 

(l-a;)(l-2/)(l-2;)' 

Proof. To obtain this density and according to Proposition [T^ it suffices to integrate the measure 
da' d(3' for a fixed x,y. Given {x,y) G the set of admissible {a',/3') G satisfies 

f (1 - x)a' + (1 - y)(3' + 1 > 0, 
or equivalently < (1 — xja' + 1 > (1 + yJ/S', 

[ (1 - y)/3' + !>(! + x)a'. 


{ xa' + y/3' + 7 = 
o' + /S' + 7 > 0, 
a' + 7 > 13', 

/S' + 7 > a', 


This domain is the interior of a triangle of vertices 0^, (^0, 7 ^), y^) shown below. 



Lemma 17. Let a,b,c E M. The area of a triangle of vertices (a, 0), (0, b) and (c, c) is ^\ab — be — 
ac\ = l\abc\\l-l-l . 

This lemma is proved by an easy determinant computation; it shows that the triangle we 
consider has area 

2 (7 + !/)(x - Vfy -T)) " (x + !/)(l - x)(l - y) ■ 

□ 

Remark 18. This density is not bounded: it tends to infinity at the three extreme points of the 
simplex. The reason is that these points are indifferent fixed points for the algorithm. This is a 
common feature for an additive continued fraction algorithm. However, and oppositely to what 
happens in dimension 1, the total mass of that density is bounded; its value is 

h f'~‘ _i_ dydx = /“ -LJSiLHdx = - 

Jo Jo (I + J/)(1-J-)(1-J/) ® Jo (x + l)(x-l) 4 
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5 The Cassaigne algorithm 

This algorithm was suggested by Julien Cassaigne (DynASS meeting, LIAFA, Paris, October 12th, 
2015)[^as a way to generate words of low factor complexity {p{n) = 2n + 1) in with arbitrary letter 
frequencies on a three-letter alphabet. An interesting aspect of this algorithm is that there are 
always only two branches. 

More formally, we can use the notations of the previous section, and define Cassaigne algorithm 
on the following partition of A = {(xi, 0 : 2 , X 3 ) G | xi > 0, X 2 > 0, X 3 > 0} up to a set of measure 
zero: 


Aa = {{Xi,X 2 ,X 3 ) G A I Xi > X3}, 
Ab = {{XI,X2,X3) G A I X3 > Xi}. 


We define two matrices: 



/ 1 1 0 \ 


/ 0 1 0 \ 

Ca = 

0 0 

0 

0 

and Cb = 

0 

0 

0 


We define the matrix function M : A —)■ GL{3,Q) such that M(x) = Ai if and only x G A* 
for i G {a,b}. Recall that this matrix function defines the functions F{x.) = M(:^“^x, /(x) = 
F(x)/||F(x)||i and F(x,a) = (M(x)“^x,M(x)'''a). We show in the left of Figure I the trace on 
the unit simplex of the partition of A; each branch of the map F sends Aj to all of A for i G {a, b}. 




Figure 6: Left: the trace on the unit simplex of the partition of A into and A;,. Right: the partition 
of A* into A* and A^. 


Explicitly, the function F is given by 


F : {x,y,z,a,f3,'y) ^ 


{x - z,z,y,a,a + -f,P) 
{y,x,z-x,/3,a + -f,-f) 


if {x,y,z) G Ai, 
if {x,y,z) G A 2 . 


^http;//www.liafa.univ-paris-diderot.fr/dyna 3 s/ 0 ct 2015 
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One can easily find a cone A* such that F is a bijection on F = A x A*: numerical experimen¬ 
tations show that after few iterations of F, {a,(3,'y) belongs to the subset of A* C A 

A* = {(a, /3, 7 ) G I max{a, 7 } < /? < a -f 7 }. 

We define a partition (see the right of Figure]^ of A* by 

= {(«>/3,7) e A* I a < 7 }, 

K = {(a,/^,7) e A* I a > 7 }. 

Lemma 19. {A*,A^} is a partition of A*, up to a set of measure 0. 

Proof. Trivial. □ 

As an immediate consequence, we obtain a partition of D: 

Corollary 20. Up to a set of measure 0, we ean write the set D = A x A* as a disjoint union 

D= Ai X A* = A X A* 

i£{a,b} 


Hence, to prove that F is a bijection on D, since each branch of F is a non-degenerate linear 
map, it suffices to prove: 

Lemma 21. We have F{Aa x A*) = A x A* and F{Ab x A*) = A x A^. 


Proof. We check directly by computation that F(Aj) = A. We want to show that if x G Aj and 
a G A*, then Cja G A*. It is enough to show that the matrix Cj sends the extreme points of A* 
to the extreme points of A* for i G {a, b}. This is done by the following matrix computation 



and this proves the equalities. See Figure 


and Cl 


/ 1 

1 


/ 1 

1 

i\ 

1 

1 

' = 

1 

2 

1 

VO 

1 

c 

lo 

1 

c 


□ 




Figure 7: The function F : Aj x A* —)■ A x A* is one-to-one and onto for each i G {1,2,3,4}. 


We now apply the results of the previous section. Recall that 

Di = {(a:, y, z, a, (3,'y) G F \ xa + yfd + Z'y = 1}. 
and the surface of section is defined by 

S = {(x, I/, z, a, /5, 7 ) G \xa + y(3 + zj = l,x + y + z = l,a<P,'y<P,P<a + 7 }. 
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It is convenient to take the coordinates on Di defined in the previous section: We keep variables x 
and y and change the coordinates a' = a — 7 , /S' = /S — 7 , e = xa + y(3 + Z''), r = — \og{x+y + z)] it is 
readily checked that the Jacobian is one. Since the domain Di is defined by e = 0, and coordinate 
r is the return time of the flow to S, we get 

E = {(x, I/, a\ I3\ T, e) G | r = 0, e = 1, a' < /S', 0 < 13', 13' < a' + 7 } 

Furthermore, the invariant measure for the return map / of the flow to E is dx dy da' dl3'. From 
this, we obtain: 

Proposition 22. The density function of the invariant measure 0 / / : A —> A for the Cassaigne 
algorithm is 

1 

(1 -x)(l 

Proof. To obtain this density and according to Proposition it suffices to integrate the measure 
da'd(3' for a fixed x,j/. Given (x, ?/) G the set of admissible {a', (3') G satisfies 


( xa' + y(3' + 'j = 1, 
I a' < /S', 

0</S', 

[ /?' < a' + 7, 


or equivalently 


r a' < /S', 

0</S', 

[ (1 - x)a' + !>(! + y)l3'. 


This domain is the interior of a triangle of vertices 


(AW) 


shown below. 



1 1 


+ y X + y 


The area of the triangle is 


1 


0 


2V X — ij \x + yJ 2(1 — x)(x + y) 

Remark 23. The total mass of that density is hounded; its value is 

1 


r*! pl — X 


lo Jo 2(1 — x)(x + y) 


dydx = 


1 log (x) , TT^ 

dx = — 
12 


'0 2 (x — 1) 


□ 


6 The Brun algorithm 

The Brun algorithm is arguably the simplest multidimensional continued fraction algorithm; it is 
defined on the positive cone A = {(xi,X 2 ,X 3 ) G | xi > 0,X2 > 0,X3 > 0} by taking the second 
largest coordinate from the largest one. Formally, we have 6 cases, and consider the following 
partition (up to a set of measure zero): 

K = {(Xi,X2,X3) G A I X^(i) < X^(2) < a:vr(3)}, 
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for 71 E S 3 where S 3 = {123,132,231,213,312,321} is the set of permutations of {1,2,3}. We 
dehne the six elementary matrices: 


-B 123 — 


-B 213 — 


1 0 0 \ 

0 10 ,5i32 = 

0 1 1 y 


1 

0 

0 \ , 

f 1 

0 

1 \ 

0 

1 

1 5 -B 23 I = 

0 

1 

0 

0 

0 

1 y ' 

1 0 

0 

1 / 


1 

0 

0 ^ 

\ 

/1 

0 

0 \ 

0 

1 

0 

5 -B 3 I 2 = 

1 

1 

0 , B 321 

1 

0 

1 y 

/ 

lo 

0 

1 ; 


1 1 0 \ 

010 

001 ; 


and the matrix function M : A —)■ GL{3, Q) is dehned according to the above partition: 


M : X I— )■ if and only if x G Aj^. 

As in Section the matrix function M (x) allows to dehne the functions 

F(x) = M(x)~^x, /(x) = and F(x,a) = (M(x)“^x, M(x)^a). 

11-^ (xjlli 

The six branches of the map F are not full; that is, their image is not the complete cone. It 
will be useful to dehne the subcones 0 ^ by 


©TT = {{Xi,X2,X3) e A I X^(1) < X^(2)}. 

Remark 24. Each is a union of three A„, for all permutation a such that cT“^ 7 r(l) < cT“^ 7 r( 2 ); 
a direct computation shows that F{At^) = 0,^, so that the A^ form a Markov partition for F, which 
can he proved to he generating. 

We want to apply the techniques of Section however, as experiments show, the domain for 
the natural extension can not be written as a global product, but as a disjoint union of products. 
Dehne for each permutation 71 E S 3 , 


^*TT — {(0^150^2,0:3) ^ A I 0 : 7 r(l) < 0 : 7 r( 3 ) < a^( 2 )} 

= {( 0 : 1 , 02 , 0 : 3 ) G A I a; 7 r(l) < Q: 7 r( 3 )} 



Figure 8 : The function F : A^^ x 0* —)■ 0,^ x A* is one-to-one and onto for each tt E S 3 . 

This is just a convenient notation; if we denote by p = (123) the transposition exchanging 2 
and 3, we have A* = A^^p and 0* = Qnp, and the previous remark applies: each 0* is the disjoint 
union of three sets A*. We hrst prove a lemma. 
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Lemma 25. We have: 


U (A. X e;) = u (0 . X a;) 

TrS^a ttScSs 

up to a set of measure 0 , where both unions are disjoint. 

Proof. Since all are disjoint, the first union must be disjoint, and since all A* are disjoint, the 
second union is disjoint. 

By the previous remark, each 0* is the disjoint union of three sets A^, for those a E S 3 such 
that cT“^ 7 r(l) < ct“V( 3). This equality is possible if, and only if, (t( 3) = vr(3), or (t(3) = 7 r( 2 ) and 
ct(2) = 7 r( 3 ). Hence the first set is the union of all the 18 products x Ayr, for all pair (cr, vr) 
satisfying this condition. 

The second set can be written Utts^s (®a- x A^) = UTrs^g (©tt x A^p) = UTrs^a (©^p ^ ^a-), where 
p = (132); and A^- C ©Trp if and only if a~^Tip{l) < (j“^ 7 rp( 2 ): by dehnition of p, this is the same 
condition, hence the two sets are equal. □ 

Proposition 26. The function F restricted to the domain A^^ x 0* is one-to-one and onto 0,^ x A* 
for each permutation n E S 3 . The function F on the domain D = U.,re 53 (^tt x 0*) is one-to-one 
and onto (See Figure^. 

Proof. We first prove that F (A^r x 0 *) = 0 ^ x A*; by symmetry, it is enough to prove it for tt = 
( 1 , 2 , 3 ). Let {xi,X2,X3,ai,a2,a3) E A^x 0 ;, and let (pi, ps, 2/3,/^i,/ 32 ,/^a) = P(xi, Xa, 0:3, «!, ^2, as) = 
(xi, X2, X3—X2, oi, oa+aa, 03); we check that yi < 1/2 and (di < (^3 < (^2, which is the condition to be 
in 0 ^ X A*. To prove equality, it is enough to consider the reciprocal map ^ (j/i, pa, 2/3, /Si, /d2, /Sa) = 
(2/I5 2/25 2/2 + 2/35 /Si, /S2 — /S3, /S3), and to prove that F~^ ( 0 ^ x A*) = A^ x 0 *, which is done exactly 
the same way. 

Since each branch of F is an invertible linear map, it is a bijection between its domain and its 
image; by the previous lemma, F is a bijection on D. 

□ 

Proposition 27. The density function of the invariant measure 0 / / : A —A for the Brun 
algorithm is 

1 

2 ^ 7 r( 2 )(f ^ 7 r( 2 ))(f ^ 7 r(l) ^ 7 r( 2 )) 

on the part x = (xi, Xa, 0 : 3 ) G A.^. fl A. 

Proof. In this proof, we suppose that vr = (1,2,3). For Brun algorithm, the surface of section for 
the part x G A^^ is 

3 3 

S = {(xi,Xa,X3,ai,aa,a3) G M+ | '^XiUt = 1, ^ x* = 1, x G A^, a G 0*}. 

^=1 i=l 

As explained in section we use an explicit change of variables that is useful to compute explic¬ 
itly all these measures in the case of the norm ||x|| = I]i=i 13^*I- Define yi = xi, 1/2 = xa and 
T = — log(X)f=i 1 ^* 1 ): define (3i = ai — 0 / 3 , /Sa = aa — 0/3 and e = A straightforward 

computation shows that this change of coordinates has jacobian 1. Below, we keep variables Xi 
and Xa instead of yi and i/a- We get that S is the set of (xi, Xa, /?i, /Sa, e) G satisfying 

r = 0, A < 0, 

e = 1 , /?! -[- 0/3 > 0 , 

xi < Xa, /S 2 + a3 > 0, 

0/3 > 0 . 
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The equalities e = 1 and ^ allows to substitute as above. Indeed, 

3 

1 = e = Q!jXj = {^l3i-\-Ol^)Xi-{-{^l32-{-Ol^)x2-\-Ol^X2, = l3iXi-\-l32X2-\-Oiz{,^l~^^2^X^) = l3iXi-\-l32X2-\-Oi‘i- 

i=l 

So we get that S is described by 


r = 0 , A < 0 , 

e = 1 , jSiXi + (52X2 - /5i < 1 , 

Xi < X2, (5iXi + (52X2 - /32 < 1 , 

(5iXi + (52X2 < 1 . 


From Proposition 11, the density is equal to the volume of the polytope {a : (x, a) G S}: 


6 {Xi,X 2 ) = / 


{(/3i,/32):(a;i,a:2,0,/3i,/32,l)eS} 


1 d(5i d(52 


Given (xi, X2) G the polytope corresponds to the set of (/Si, /S 2 ) ^ 1^^ satisfying the four inequal¬ 
ities. This domain is the interior of a triangle with vertices (^0, (o, 

shown below. 



Its area is 


2 \X2 X 2 -I 


0 - 


a^i+a; 2 -l/ 2x2(1 - 0:2 )(1 - xi - 0:2) 


□ 


The mass of that density for the part Xi < X2 < X 3 is: 


,1 „ 1 XI 


l(J( ' 2x,(l - x,)(l - X. - = -^+'2 


we will prove in Section that this quantity is equal to 


TT^ 

24* 
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Fl] 


• Xn+'l 


/ •Xn 


Figure 9: For some starting point (x, a) G xM*^, we consider the sequence (x„, a„) = F"”(x, a) for n > 0. 
The points x„, a„, x„_|_i and a„+i are drawn in four different plots with a color or label vr associated to 
the matrix M(x„). This experimentation allows to view the domain of the natural extension of F. 
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Figure 10: Experimental natural extension of Brun algorithm. 

7 Experimentations: domain of the natural extension 

Since the map F is linear, the trajectory of Ax is homothetic to that of x; hence, we can always 
normalize x to and consider it as a point in the nnit simplex. We can do the same thing 

for a, and both (x,a) G x and F(x,a), its image nnder F, can be represented graphically 
by a pair of points in a pair of nnit simplexes. This gives fonr points in four distinct simplexes. 
To visualise the effect of F, we draw each of these four points with the color according to which 
element of the partition of A the point x belongs. We could do this exercice for a set of random 
points (x, a) G X but in order to hnd the domain of the natural extension, it is better to 
consider some orbit (F"'(x, a))„>o for some starting point (x, a) G x M'^. Below are the pictures 
we get for some Multidimensional Continued Fraction algorithms. 

The experimentation shows that a domain of positive measure can be expected for Brun and 
Reverse algorithms. The method also works for Selmer algorithm. The method does not work for 
Poincare and Meester algorithms: both do not have a hnite measure preserving transformation. 



Figure 11: Experimental natural extension of Reverse algorithm. 
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Figure 12: Experimental natural extension of Cassaigne algorithm. 



Figure 13: Experimental natural extension of Arnoux-Rauzy-Poincare algorithm. 


Note that Poincare does have an invariant measure equivalent to Lebesgue but it is not conservative. 
Finally, the domain for AR-Poincare BL15 and for the additive version of Jacobi-Perron algorithm 
are fractals and dehned by an IFS. Both of these algorithms are known to have an invariant 
measure absolutely continuous with respect to Lebesgue measure; this was shown for Jacobi- 
Perron algorithm by Broise and Guivarc’h BAGOl , and we will prove it for AR-Poincare in a 


forthcoming paper, but we do not know their density function. 

In Figure we compute the hrst 2000000 iterations of an orbit in a 1024 x 1024 pixels picture. 

This set presents a symmetry of order 3; it is decomposed in 9 subsets, and one can prove that 
there are 3 disjoint groups of three sets. It is unclear whether these sets contain open balls, have 
positive measure, or even dimension 2; it is also unclear if the three sets in each groups are disjoint, 
but if this is the case, it seems difficult to prove. 


8 Additional remarks 


8.1 Invariant domain and fixed point theorem 

Of course, the essential element for this construction is to hnd the invariant domain D. It can 
sometimes be proved from elementary considerations, as we have seen in the two examples above, 
and numerical experimentations can be useful. 

When no obvious solution is found, one can use the fact that the formula for / is known. For 
hxed X, / is an affine map on an hyperplane in A*. In some cases, it can be proved that this 
collection of maps depending on x is uniformly cont ractin g, and this implies that there exists a 
unique compact invariant set which is hxed by /, see AS15 . Numerical experiments seem to show 


that this property is still valid under much weaker conditions. This allows, as we have seen above, 
to guess the domain. 

However, this domain is interesting only if it has strictly positive measure, which is not clear 
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Figure 14: 2 x 10® iterations shown in a 1024 x 1024 pictures of the window [—0.6,0.6] x [—0.6, 0.6]. On the 
left, points in the six Poincare branches are drawn last. On the right, points in the three Arnoux-Rauzy 
branches are drawn last. 


in a number of cases; even in that case, it can be very difficult to describe, as shown by the case of 
the Japanese continued fraction, which has been very precisely studied, see for example 


KSS12 


Remark also that a very well studied case, that of Rauzy induction, show the complications 
which can appear : if one uses exactly the procedure we describe in the present paper, one does 
not get a natural extension, because one gets a subset of lower dimension (given by the subspace 
called H{ti) by Veech, see Vee82| ), and it is necessary to add other coordinates to get the natural 
extension. 


8.2 Symmetric algorithms 

The group Sd of permutations on d elements acts in a natural way on by <J.{xi,... ,Xd) = 
{xa{i), ■ ■ ■ ,Xa{d))- We will say that a map F is symmetric if, for all permutations a, F{a.x.) = 
a.F{x.). Brun and Reverse algorithms are symmetric; a classical example of a non-symmetric 
algorithm is Jacobi-Perron algorithm: a permutation of coordinates gives in general a very different 
orbit. 

When we have a symmetric algorithm dehned on the positive cone, we can quotient by Sd, and 
restrict the algorithm to the subset dehned hy 0 < xi < X 2 <■■■< Xd- The original symmetric 
algorithm is often named the unsorted algorithm, and the quotient algorithm is named the sorted 
algorithm. These algorithms are essentially equivalent, since the sorted algorithm is a dl-to-l factor 
of the unsorted one. The sorted algorithm is generally more convenient for explicit computations, 
in particular numerical experiments, since the relations between coordinates are known; but the 
unsorted algorithm, being more symmetric, gives usually simpler and more natural formulas, and 
is better adapted to theoretical considerations; another interesting feature is that, in classical 
examples, the unsorted algorithm preserves the orientation, but the sorted algorithm does not. 
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Figure 15: Left:10® iterations and zoom in [0.05,0.15] x [0.05,0.15]. RightilO® iterations and zoom in 
[0.09,0.11] X [0.09,0.11]. 


since transpositions reverse the orientation. 

This leads to many apparently distinct, but essentially equivalent, presentations for the same 
algorithm; this is already true for the basic Farey algorithm, which takes several forms with similar 
properties. 

8.3 Parabolic points and acceleration 

We say that a continued fraction algorithm has an indifferent fixed point if one of the branches 
of the projective map / associated to this algorithm has a hxed point where the Jacobian matrix 
is the identity (this does not depend on the choice of coordinates). Several classical algorithms 
have one or several indifferent hxed points; this presents the inconvenient, in dimension 1, that 
the invariant measure becomes inhnite, and the entropy zero. 

A well-known way to get rid of this problem consists in "accelerating" the algorithm; more 
precisely, let fi be the branch of / which has an indifferent hxed point, and let Ai be its domain. 
For any x diherent of the hxed point, dehne Ux as the smallest n such the f^{x) ^ A^, and replace 
the branch fi by the countable collection of branches x h-)■ f^^{x). This replaces what is generally 
called an "additive" algorithm with a hnite number of branches by a "multiplicative" algorithm, 
with a countable number of branches, no indiherent hxed point, and with a hnite invariant measure. 

8.4 Changes of coordinates 

As we said before, a continued fraction map in the usual sense is a projective map; hence there is 
no canonical coordinates for explicit computations, and many things depend in particular on the 
choice of the section E. As an example, we show another section for Brun’s algorithm, which will 
allow us to compute the total mass of the invariant measure in a diherent way. 
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Instead of using the norm I]i=i ..3 chose the norm supj^^ 3 |xi|, and we work out the 

proof of Proposition [2^ in this setting. We suppose as before that tt = (1, 2, 3). The corresponding 
section, for this choice of norm, is 


3 

S = {(xi, X 2 , X 3 , «!, 02, tt3) e M+ I ^ XiOi = l,Xi < X 2 < X 3 = 1, «! < Os}. 

i=l 

We make a different choice of coordinates here, adapted to this norm : keep xi, X 2 , ai and 02 ; 
define r = — log(|a: 3 |) and e = A straightforward computation shows that this change 

of coordinates has Jacobian 1 . 

We get that S is the set of (xi, X 2 , (^ 2 , r, e) G satisfying 


j 0 < xi < a ;2 < 1 , r = 0 , e = 1 

I 0 < «! < 03, 0 < a2- • 

The equalities e = 1 and X 3 = 1 allows to substitute above, hence S is also described by 


0 < Xi < a ;2 < 1 , r = 0 , e = 1 

cri(l T Xi) T 02X2 ^ 1 ) 0 < CTi, 0 < 02- 


Hence the density is equal to the volume of the polytope {a : (x, a) G S}: 


( 5 (Xi,X 2 ) = / 

^ I 


{(ai,a2)-.(xi,X2,0,oii,a2,l)&'^} 


1 daida2 


This is again given by a triangle, this time with vertices (0, 0), (1/(1 + Xi), 0), (0, 1 /X 2 ). Hence the 
invariant density is given in these coordinates by 2 x 2 (i+xi) ~ 

The mass of that density on the triangle 0 < xi < X 2 < 1 is: 



1 

2 x 2(1 + Xi) 


dxidx2 


log(l + X 2 ) 

Jo 2X2 


Taking the series expansion of log(l + X 2 ) and integrating, we find that this quantity is equal 

/ 1 \ 71 1 2 

to the series " 2^2 = fj, as announced at the end of section 5. This is just a change of 

variable in the integral, which can be made completely explicit. Since we made the computation 

2 

on a sixth of the surface, the total mass of the invariant measure for the unsorted algorithm is 
The equality —+ 1 log (3) log (2) — 1 Li 2 (|) + | Li 2 (|) — Li 2 (—|) = |^ is probably well-known 
among the many identities for the dilogarithm function, but we could not find it in the literature. 


8.5 Variants of the classical continued fraction 

We have studied in Section the unsorted Farey map. One of the advantages of this map is that it 
preserves orientation, and is related to matrices in SL{2,'R). However, most presentations prefer 
to deal with the sorted map, restricting to the cone x < y, and in that case it is convenient to 
take as section y = 1 ; this leads to the usual sorted Farey map, given by x h->■ if x < 1 and 

xi-G-l — lifx>l. 

X 2 

This map has only one indifferent fixed point, at 0. Taking the acceleration of the map at this 
fixed point leads to the usual Gauss map, x h-)■ 

Many other variants are possible, by taking different sections; they are all related to the geodesic 
flow on the modular surface and the group SL{2,Z). 
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8.6 Higher dimensions 

It is natural to generalize to higher dimensions. For the reverse algorithm, a direct generalization 
is not possible: the basis of the algorithm is that the complement in the unit simplex of the three 
simplexes of size | given by the condition Xj > | is also a simplex, which can be sent to the whole 
simplex by an homothety of ratio —2. In dimension 4, the complement, in the unit simplex, of the 
four simplexes x* > | is no more a simplex, but an octahedron, so we cannot extend our simple 
solution to that case. 

But it is possible, and easy, to generalize Brun algorithm to any dimension. We dehne it by 
subtracting the second largest coordinate by the largest coordinate. We can dehne an invariant 
domain for the natural extension in the same way. Let A be the positive cone in M'^, and dehne, 
for tie Sd 


= {(xi,..., Xrf) e A I x^(i) < x^( 2 ) < ■ ■ ■ < x^(rf)} 

Ott . . . , ^ ^ I ^7r(l) ^ ^7r(2) ^ ^ ^7r(d—1)} 

{(^1; • • • : ^d) ^ A I ^n(d) ^ 1) If ^ ^ d 1} 

©tt {(®1) ■ • • 1 ®(i) ^ A I ttTrp) < Oiitid) if ^ ^ d Ij’ 


The same proof as in dimension 3, with heavier notations, proves that F sends A^ x 0; to 
0^ X A*. It is then straightforward to compute the invariant measure; similar formulas were given 
in 


AN93 and HKOl 


Proposition 28. The density function of the invariant measure 0/ / : A —^ A for the d- 
dimensional Brun algorithm is 


d-l 


(5(x) = 


(d !)• ^ J2i£Ak^i 


on the part x G A.,r where Ai = {7r(d — 1)}, Ad-i = {1, 2,..., d} \ {7r(d)} and \Ak\ = k. 
For example, taking tt = Id, when d = 2: 


(5(xi) = 


Xi (1 - Xi) 


When d = 3: 

When d = 4: 


d(xi,X2) = 


2x2 (1 - xi - a:2)(l - X 2 ) 


1 


d(Xi,X2,X3) = 


+ 


1 


(1—a:i—X3)(l— 2 : 3 ) (1—X2—X3){1—X3) 

6 X3(1 -Xi-X2 -X3) 


When d = 5, d(xi, X2, X3, X4) is 

_1_L_ 1 _ _ 1 _I_ 1 _ _ 1 _I_ 1 _ 

(1-X2-X4)(1-X4) , (l-rcj-a:4)(l-a:4) (1-3:3-3:4) (1-^4) , (I-X2-0^4 )(1-^^4) (1-X3-3:4) (1-3:4) 

I — X1—X2—X4 ' 1 — Xi—X^—X4 ' I—X2—XS—X4 

24x4 (1 — 3^1 — X2 — X3 — X4) 
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Proof. In this proof, we suppose that tt = [1,2,... ,d). For Brun algorithm, the surface of section 
for the part x G A,r is 

d d 

S = {{xi,... ,Xd,ai,... ,ad) e M+ | XI = l,'^Xi = l,:>c e A^,a e 0*}. 

i=l i=l 

As explained in section]^ we use an explicit change of variables that is useful to compute explicitly 
all these measures in the case of the norm ||x|| = Yfi=i \xi\. Dehne Hi = Xi for i < d, and yd = 
log(I]f=i k*l)j dehne (3i = ai — ad ioi i < d and /3d = J2i=iaiXi. A straightforward computation 
shows that this change of coordinates has Jacobian 1. Below, we keep variables Xi for i < d instead 
of yi. We get that E is the set of (xi,..., Xd-i, yd, Pi, ■ ■ ■, Pd-i, Pd) G satisfying 

f 2/d = 0, /Si < 0 for i < d - 1, 

I Pd = I, /Si + Qfrf > 0 for z < d, 

[ xi< ■■■ < Xd-i, ad > 0. 


The equalities /S^ = 1 and I]f=i Xi = 1 allows to substitute ad above. Indeed, 


d d— 1 d— 1 d d— 1 

^ = Pd = = ^{Pi + Old)Xi + adXd = X Pi^i + ttrf X H Pi^i + ^d- 

i=l i=l 2=1 2=1 2 = 1 


So we get that S is described by 


2/d = 0, A < 0 for z < d - 1, 

/3d = 1, Ej=i PjXj - Pi < I ioT i < d, 

Xi < • • • < Xd-l, E?=l PjXj < 1. 


From Proposition 11, the density (5(x) is equal to the volume of the polytope {b : (x, b) G S}: 


d(x) = 


'{b:(x,b)eS} 


1 db 


Given x = (xi,..., Xd-i) G \ the polytope corresponds to the set of b = (/Si,..., Pd-i) G ^ 
satisfying the above inequalities. We conclude the proof by using Lemma □ 


Lemma 29. Let I C {1, 2, ..., d} and i ^ I. Let (xi, X 2 ,..., Xd) G be positive real numbers 
such that YPi=iXi = 1. Let be the \I\-dimensional polytope defined by the set of real vectors 
{Pi : i E /) that satisfy: 

( Pi <0 fori e /\ {^}, 

< Y/jei PjXj - Pi <l fori e I, 

[ J2jei PjXj < 1. 


Then 


Vol{Pd^,) 


1 

im-Ei^i^i) 


Y. 


Proof. First notice that Pi = , ^ for all z G / is a solution to the second equation. This is a 

vertex of the polytope that we denote S. The interior of the polytope can be partionned into |/| — 1 
part pk according to the index k E I \ {£} for which P^ > Pi for all z G /\ {£}. Note that the vertex 
S is also a vertex of each of these parts. In fact, we can see that each part is in fact a pyramid 
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with top vertex S and base given by the equation 13k = 0. That base is a (|/| — l)-diniensional 
polytope given by the equations 


( l3i < 0 ioT i E I \ {k, i}, 

I EjeAtfc} - Pi < 1 for i E I\ {k}, 

{ J2j&i\{k} Pj^j < 1 

which is exactly Pi\{k},e and whose volume is Vol{Pi\{k},e)- □ 
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